
// general math test routine

#include <vector>
#include "Math.hpp"
#include "Vector.hpp"
#include "Matrix.hpp"
#include "Poly.hpp"

//GNU multiple precision
#include <gmpxx.h>

#include <iostream>

using namespace std;

int
main (int argc, char** argv)
{
  double V0[] = { 1.1, 2.0, -3, 4 };
  double V1[] = { 0,   1,    4, 2 };
  double x;
  float y;
  int j;

  DPoly P0(V0,3);
  DPoly P1(V1,3);

  if (argc > 1)
    x = atof(argv[1]);
  else
    x = 0;
  y = x;

  cout << "P0" << endl;
  cout << P0;
  cout << "P1" << endl;
  cout << P1;

  DPoly Pd = P0.Deriv();
  DPoly Pi = P0.Integral(0.0);
  cout << "P0 derivative" << endl;
  cout << Pd;
  cout << "P0 integral" << endl;
  cout << Pi;
  cout << "Pi(" << x << ") = " << Pi(x) << endl;

  DPoly P2 = P0 * P1;
  cout << "P2" << endl;
  cout << P2;

  cout << "P0(" << x << ") = " << P0(x) << endl;
  cout << "P1(" << x << ") = " << P1(x) << endl;
  cout << "P2(" << x << ") = " << P2(x) << endl;

  cout << endl;

  try {
    Legendre<double,double> *LP = new Legendre<double,double>[6];
    Vector<double> *VR = new Vector<double> [6];
    cout.precision(6);
    cout << "Legendre Coefficients" << endl;
    for (j=0; j < 6; ++j) {
      LP[j] = Legendre<double,double>(j);
      cout << "order " << j << ": " << LP[j];
    }
    for (x = -1; x <= 1; x+=1) {
      for (j=0; j < 6; ++j)
	cout << "LP" << j << "(" << x << ") = " << LP[j](x) << endl;
      cout << endl;
    }
    for (j=1; j < 6; ++j) {
      cout << "Roots for order " << j << endl;
      VR[j] = LP[j].Roots();
      cout << VR[j];
      cout << "value at roots" << endl;
      for ( int k=0; k < VR[j].size(); ++k )
	cout << LP[j](VR[j][k]) << " ";
      cout << endl;
      
    }

    // multiple precision polynomial
    //vector<mpf_class> VMP(4) = { 1.1, 2.0, -3, 4 };
    Legendre <mpf_class,mpf_class> MPF(5);

    cout << "Multiple Precision" << endl;
    cout << "Coefficients" << endl;
    for (int i=0; i <= 5; ++i)
      cout << MPF.Coeff(i) << " ";
    cout << endl;

    cout << "finding roots" << endl;
    Vector<mpf_class> roots = MPF.Roots();

    for ( int i=0; i < roots.size(); ++i )
      {
	cout << roots[i] << ", " << MPF(roots[i]) << endl;
      }
    cout << endl;

    Tchebycheff <double,double> *TP = new Tchebycheff<double,double>[6];
    cout << "Tchebycheff Coefficients" << endl;
    for (j=0; j < 6; ++j) {
      TP[j] = Tchebycheff <double,double>(j);
      cout << "order " << j << ": " << TP[j];
    }
    for (j=1; j < 6; ++j) {
      cout << "Roots for order " << j << endl;
      VR[j] = TP[j].Roots();
      cout << VR[j];
      cout << "value at roots" << endl;
      for ( int k=0; k < VR[j].size(); ++k )
	cout << TP[j](VR[j][k]) << " ";
      cout << endl;
    }
    cout << endl;

    Hermite <mpf_class,mpf_class> *HP = new Hermite<mpf_class,mpf_class>[6];
    cout << "Hermite Coefficients" << endl;
    for (j=0; j < 6; ++j) {
      HP[j] = Hermite <mpf_class,mpf_class>(j);
      cout << "order " << j << ": " << HP[j];
    }
    for (j=1; j < 6; ++j) {
      cout << "Roots for order " << j << endl;
      Vector<mpf_class> HR = HP[j].Roots();
      cout << HR;
      cout << "value at roots" << endl;
      for ( int k=0; k < HR.size(); ++k )
	cout << HP[j](HR[k]) << " ";
      cout << endl;
    }


  }
  catch (exception& ex) {
    cerr << ex.what() << endl;
  }
  catch (...) {
    cerr << "unexpected exception" << endl;
  }
}

